###############################################################
## The following code is part of the replication archive of  ## 
## Genovese, Kern & Martin (2016) "Policy Alternation", ISQ  ## 
###############################################################

#install.packages("systemfit")
#install.packages("stargazer")
#install.packages("nls2")
#install.packages("visreg")
library(nlme)
library(nls2)
library(lattice)
library(systemfit) 
library(foreign)
library(stargazer)
library(coefplot)
library(effects)
library(colorspace)
library(splines)
library(ggplot2)
library(visreg)
library(systemfit) 
#library(reshape2)
#library(readstata13)

# Read the data into R:
scr <- read.dta("~/Dropbox/.../gkm_replication_isq/carbon policy/EU neighbors_carbon policy.dta")
scr$countryn<-as.factor(as.character(scr$countryn))

scr$l_sl_co2tax<-lag(scr$sl_co2tax)*100
scr$l_sl_co2allowances <-lag(scr$sl_co2allowances) #already adjusted

##################
# OLS: Table A.6 #
##################

co2tax     <- co2taxnormal  ~   lag_co2taxnormal + allowances_gdp_lag + gdppercapita + gdppercapitasq + co2ktpercapita + co2ktpercapitasq + energydependencektofoileq + govteffectiveness + euintegration + dpi + ratio_over_eutob + sl_co2tax   
co2trade     <- allowances_gdp ~   allowances_gdp_lag + lag_co2taxnormal + gdppercapita + gdppercapitasq + co2ktpercapita + co2ktpercapitasq + energydependencektofoileq + govteffectiveness + euintegration + dpi + ratio_over_eutob + sl_co2allowances 
labels<-list("co2tax", "co2trade") 
system<-list(co2tax, co2trade) 
fitolsa <- systemfit( system, data=scr )
summary(fitolsa)
out <- capture.output(summary(fitolsa))

co2tax     <- co2taxnormal  ~   lag_co2taxnormal + allowances_gdp_lag + gdppercapita + gdppercapitasq + co2ktpercapita + co2ktpercapitasq + energydependencektofoileq + govteffectiveness + euintegration + dpi + ratio_over_eutob + sl_co2tax  + l_sl_co2allowances 
co2trade     <- allowances_gdp ~   allowances_gdp_lag + lag_co2taxnormal + gdppercapita + gdppercapitasq + co2ktpercapita + co2ktpercapitasq + energydependencektofoileq + govteffectiveness + euintegration + dpi + ratio_over_eutob + sl_co2allowances + l_sl_co2tax
labels<-list("co2tax", "co2trade") 
system<-list(co2tax, co2trade) 
system<-list(co2tax, co2trade) 
fitolsb <- systemfit( system, data=scr )
summary(fitolsb)
out <- capture.output(summary(fitolsb))

co2tax   <- co2taxnormal  ~ lag_co2taxnormal + allowances_gdp_lag + gdppercapita + gdppercapitasq + co2ktpercapita + co2ktpercapitasq + energydependencektofoileq + govteffectiveness + euintegration + dpi + ratio_over_eutob + sl_co2tax  + l_sl_co2allowances + I( l_sl_co2allowances* ratio_over_eutob)    # 
co2trade  <- allowances_gdp  ~ allowances_gdp_lag + lag_co2taxnormal + gdppercapita + gdppercapitasq + co2ktpercapita + co2ktpercapitasq + energydependencektofoileq + govteffectiveness + euintegration + dpi + ratio_over_eutob + sl_co2allowances + l_sl_co2tax+ I( l_sl_co2tax* ratio_over_eutob)  # 
labels<-list("co2tax", "co2trade") 
system<-list(co2tax, co2trade) 
fitolsc <- systemfit( system, data=scr )
summary(fitolsc)
out <- capture.output(summary(fitolsc))

###################
# 3SLS: Table A.7 #
###################

co2tax   <- co2taxnormal  ~ lag_co2taxnormal + allowances_gdp_lag + gdppercapita + gdppercapitasq + co2ktpercapita + co2ktpercapitasq + energydependencektofoileq + govteffectiveness + euintegration + dpi + ratio_over_eutob + sl_co2tax # 
co2trade  <- allowances_gdp  ~ allowances_gdp_lag + lag_co2taxnormal + gdppercapita + gdppercapitasq + co2ktpercapita + co2ktpercapitasq + energydependencektofoileq + govteffectiveness + euintegration + dpi + ratio_over_eutob + sl_co2allowances#  
labels<-list("co2tax", "co2trade") 
system<-list(co2tax, co2trade)                                                       
inst1<- ~   lag_co2taxnormal +  gdppercapita + gdppercapitasq + co2ktpercapita + co2ktpercapitasq + energydependencektofoileq + govteffectiveness + euintegration + dpi + ratio_over_eutob + sl_co2tax + countryn
inst2<- ~  allowances_gdp_lag +  gdppercapita + gdppercapitasq + co2ktpercapita + co2ktpercapitasq + energydependencektofoileq + govteffectiveness + euintegration + dpi + ratio_over_eutob  + sl_co2allowances  + countryn                                                                                
inst<-list(inst1, inst2)
fit3sls_a<-systemfit( system,"3SLS",  inst=inst,  data = scr )  
summary(fit3sls_a)  

co2tax   <- co2taxnormal  ~ lag_co2taxnormal + allowances_gdp_lag + gdppercapita + gdppercapitasq + co2ktpercapita + co2ktpercapitasq + energydependencektofoileq + govteffectiveness + euintegration + dpi + ratio_over_eutob + sl_co2tax + l_sl_co2allowances# 
co2trade  <- allowances_gdp  ~ allowances_gdp_lag + lag_co2taxnormal + gdppercapita + gdppercapitasq + co2ktpercapita + co2ktpercapitasq + energydependencektofoileq + govteffectiveness + euintegration + dpi + ratio_over_eutob + sl_co2allowances + l_sl_co2tax#  
labels<-list("co2tax", "co2trade") 
system<-list(co2tax, co2trade)                                                       
inst1<- ~   lag_co2taxnormal +   gdppercapita + gdppercapitasq + co2ktpercapita + co2ktpercapitasq + energydependencektofoileq + govteffectiveness + euintegration + dpi + ratio_over_eutob + sl_co2tax + countryn
inst2<- ~  allowances_gdp_lag +   gdppercapita + gdppercapitasq + co2ktpercapita + co2ktpercapitasq + energydependencektofoileq + govteffectiveness + euintegration + dpi + ratio_over_eutob  + sl_co2allowances  + countryn                                                                                
inst<-list(inst1, inst2)
fit3sls_b<-systemfit( system,"3SLS",  inst=inst,  data = scr )  
summary(fit3sls_b)  

co2tax   <- co2taxnormal  ~ lag_co2taxnormal + allowances_gdp_lag + gdppercapita + gdppercapitasq + co2ktpercapita + co2ktpercapitasq + energydependencektofoileq + govteffectiveness + euintegration + dpi + ratio_over_eutob + sl_co2tax  + l_sl_co2allowances + I( l_sl_co2allowances* ratio_over_eutob)    # 
co2trade  <- allowances_gdp  ~ allowances_gdp_lag + lag_co2taxnormal + gdppercapita + gdppercapitasq + co2ktpercapita + co2ktpercapitasq + energydependencektofoileq + govteffectiveness + euintegration + dpi + ratio_over_eutob + sl_co2allowances + l_sl_co2tax+ I( l_sl_co2tax* ratio_over_eutob)  # 
labels<-list("co2tax", "co2trade") 
system<-list(co2tax, co2trade)                                                       
inst1<- ~  lag_co2taxnormal +   gdppercapita + gdppercapitasq + co2ktpercapita + co2ktpercapitasq + energydependencektofoileq + govteffectiveness + euintegration + dpi + ratio_over_eutob + sl_co2tax + countryn
inst2<- ~  allowances_gdp_lag +   gdppercapita + gdppercapitasq + co2ktpercapita + co2ktpercapitasq + energydependencektofoileq + govteffectiveness + euintegration + dpi + ratio_over_eutob  + sl_co2allowances  + countryn                                                                                
inst<-list(inst1, inst2)
fit3sls_c<-systemfit( system,"3SLS",  inst=inst,  data = scr )  
summary(fit3sls_c)  

#co2tax   <- lm(co2taxnormal  ~ lag_co2taxnormal + allowances_gdp_lag + gdppercapita + gdppercapitasq + co2ktpercapita + co2ktpercapitasq + energydependencektofoileq + govteffectiveness + euintegration + dpi + ratio_over_eutob + sl_co2tax  + l_sl_co2allowances + I( l_sl_co2allowances* ratio_over_eutob) , data =scr )  # 
#visreg2d(co2tax, y="l_sl_co2allowances", x="ratio_over_eutob",plot.type="persp", zlab="CO2 tax", ylab="CO2 allowances spatial lag", xlab="actual economic flows", cex.lab=.9, cex.axis=.8, cex.main=.5, cex.sub=.6)

